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A MODIFIED DODGE ALGORITHM FOR THE PARABOLIZED NAVIER-STOKES 


EQUATIONS AND COMPRESSIBLE DUCT FLOWS 
By 

C. H. Cooke* 

SUMMARY 

A revised version of Dodge's split-velocity method for numerical 
calculation of compressible duct flow has been developed. The revision 
incorporates balancing of mass flow rates on each marching step in ordtir 
to maintain front-to-back continuity during the calculation. The (checker- 
board) zebra algorithm is applied to solution of the three-dimensional 
continuity equation in conservative form. A second-order A-stable linear 
multistep method is employed in effecting a marching solution of the para- 
bolized momentum equations. A checkerboard SOR iteration is used to solve 
the resulting implicit nonlinear systems of finite-difference equations- 
which govern stepwise transition. 


INTRODUaiON 

It has been said that the full Navier- Stokes equations represent the 
ultimate mathematical model upon which to bass numerical algorithms for 
predicting flows of practical significance. However, even with the advent 
of the so-called vector computers with vast virtual memory and quadrupled 
processing speeds, extant numerical and computational difficulties are 
sufficient to merit a search for simpler mathematical models and less 
complicated numerical methods which can still provide useful solutions to 
problems of interest. Thus, considerable analysis and numerical experiment 
has been devoted to the exploitation of parabolized marching methods for 
flow prediction. References 1 to 7 represent a perhaps typical but by 
no means exhaustive sampling of the available literature on this subject. 
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Norfolk, Virginia 23508. 


The parabolized inarching methods are somewhat more general in applica- 
tion than the classical boundary- layer approach, since transverse pressure 
gradients are not disregarded and, in some cases, upstream influences 
can be transmitted through the pressure field. However, the basic 
assumption that streamwise viscous diffusion can be neglected restricts 
application to flows with a primary flow direction, limited upstream 
influence, and which may exhibit at worst crossplane recirculation. 
Unfortunately, in subsonic and transonic wind-tunnel flows, the elliptic 
upstream influence can be a significant factor in the flow dynamics; 
hence, interest arises in simpler mathematical models which permit this 
interaction. A case in point has been the development of Dodge's velocity 
splitting method, which allows global propagation of influence through the 
pressure field, and which has met with successes in both unconfined com- 
pressible and confined incompressible flows (ref. 7-10). However, as yet 
the method is by no means fully proven. 

In this paper we shall be concerned with the application of a com- 
pressible formulation of Dodge's split velocity technique (ref. 9) to the 
calculation of developing (nonentry region) flow in a square duct. The 
original method has been revised to effect constant mass flow rate on each 
transverse plane while marching down the channel. Parabolized momentum 
equations are employed. However, a fully elliptic pressure field is allowed 
by the iterative manner in which the solution of the continuity equation 
is coupled into the calculation procedure. Application of the presently 
developed computer algorithm is restricted to subsonic flow. It could 
readily be altered to allow transonic calculations through modification 
or replacement of the algorithm used to solve the conservative continuity 
equation. 

Computational simplicity as well as numerical stability is achieved 
in marching the momentum equations with' an A-stable (ref. 11) implicit 
linear multistep method, the equations of which are iteratively solved at 
each step by employing checkerboard successive overrelaxation. While 
this solution procedure may be considered expensive, the presence of 
quadratic as well as higher order nonlinearities in the parabolized 
momentum equations requires that some iteration be employed to improve 
accuracy. As an extra benefit, the wide-ranging stability of the 


resulting inarching equations appears well worth the cost. 

Finally, the peak efficiency of the methods developed is undoubtably 
best realized on the computer system for which it has been designed, 
namely, the Cyber 203. For such machines, a numerical algorithm must 
effectively exploit the array processing capabilities; otherwise, methods 
which are not highly vectorizable misuse the available computing potential 
and can result in quite ordinary processing speeds. The explicit nature 
of the checkerboard algorithm yields a highly vectorizable method ideally 
suited for the array processor. 

In certain parabolize^ m,'».*ching schemes for confined flow (ref. 1) 
it has been the practice to aecouple streamwise and transverse pressure 
gradients. Some argue (ref. 12) that this is necessary in order to obtain 
meaningful physical solutions with parabolized equations. While results 
are still -inconclusive, computational experience gained in the current 
research appears to support this belief. Weak, but not total, uncoupling 
of the streamwise pressure gradient has appeared necessary, although this 
may stem from the manner in which local continuity of mass flow is enforced. 
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PARABOLIZED GOVERNING EQUATIONS 


The nondimens ionali zed Navier-Stokes equations for compressible 
steady flow with which we shall be concerned are now written. 


Continuity; 

V • pw a 0 
Momentum; 

“(“ • * ’(j ' • ") 

T a T^ - i w2 
O 2 

Here, for flow in ducts with nonconducting walls, the usual energy 
equation has been replaced by the algebraic constant total temperature, 
relation (3) . The constitutive relations are 


( 1 ) 

( 2 ) 

(3) 


P a ^ pT (4) 

and the viscosity approximation 

w “ (Y - 1)T. (5) 

For subsonic flow the governing equations are elliptic. However, 
a common approximation used to parabolize these equations (refs. 1,2) 
is obtained by neglecting streamwise diffusion terms in equation (2). 

With the exception of the entry region, the approximation is considered 
valid for flow in channels whose lengths are large compared to half- 
width (ref. 2). Perhaps it shou'd be remarked that, when Dodge's method 
is applied in obtaining numerical solutions of these equations, the 
approximation is only a partial parabolization since the pressure field 
is obtained from an elliptic boundary value problem. This, of course, 
allows global propagation of disturbances, through the pressure field and 
the iteration process. 


DODGE'S METHOD 
Introduction 

In Dodge's method, the total velocity vector w is arbitrarily 
decomposed as a sum of rotational and irrotational parts: 
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w ■ V(> u 


( 6 ) 


where ^ is a scalar potential. Pressure is hypothesized to depend solely 
upon the irrotational velocity acc .ding to the isentropic relationship 

P « Pq(1 - (7) 


However, density is decomposed as the sum of a viscous contribution 
0 ^ and an isentropic contribution p*: 

P » Py ♦ P* (8) 

where 

p* - Pq(1 - (9) 

Substituting equations (6), (7), and (9) in equations (1) and (2) leads 
to the so-called split equations 

V • pV^ = -7 • pu (10) 


and 


p(w • 7)w - p*(7(^ • 7)7^ 7x1 


• ’(j ’ • “) 


Equations (10) and (11) are to be iteratively solved: equation (10) 

» 

with a 3-D relaxation method for elliptic equations, following which the 
parabolized version of equation (11) is marched downstream by employing 
a checkerboard iteration to solve an implicit system of equations 
at each step. A synopsis of the iteration procedure is now presented. 


( 11 ) 


Overview of Iteration Procedure 

(1) Determine a suitable initial pressure distribution Pq by estimating 
a global distribution. In this investigation, pressure on the 
first pass is assumed to be a function of only streamwise displacement, 
and a mass-balancing operation establishes the initial pressure field. 

(2) Employing the current pressure field, march a parabolized version of 
equation (11) down the duct, simultaneously storing the right-hand 
side of equation (10). [See also eq. (17)]. 

(3) Solve equation (10) [or eq. (17)] to obtain an updated pressure field. 

(4) Repeat the computational pass consisting of steps (2) and (3) until 
sufficient passes and a converged pressure field are obtained. 
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Dodge's Method Revised 

Dodge (ref. 9) reports problems arising from adjustment of front-to- 
back continuity requirements with an iteration which is similar to that 
previously outlined. It is expected that this slow convergence stems 
from incomplete satisfaction of the continuity equation which could, 
for example, be solved after the momentum march terminates in some form 
such as 

7 * ^'^“^n-1 

Hiis is in contrast to the usual parabolized marching methods for which 
both mass and velocity variables are updated at each marching step. 

Physically, in order to maintain continuity in a channel flow, the 
mass flow rate 

must remain constant at each transverse plane. However, in Dodge's 
(unrevised) method this provision is only weakly incorporated through 
equation (10), which is solved globally upon termination of a marching 
pass. Thus, poor satisfaction of mass balancing during the momentum 
marching process is only to be expected, as numerical experimentation 
indicates. 

Consequently, we have chosen to revise the Dodge technique in a manner 
which alleviates this difficulty. This was at first attempted by employing 


♦ 


* ♦ 4>(x,y,z) 

•'O 


(14) 


to write equation .(12) in the form 

V . . -7 . ^ 

The function g is determined by iterating the numerical counter- 
part of the parabolized equation (11) at each fixed marching step until 
numerical balance of mass flow rate is achieved. This is accomplished 
through gradual changes in streamwise velocity, pressure, and density 
effected by the equation 


(IS) 
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(16) 


= g’' - o.[u, - Xrf formal) * [j^pVdz] 

with ot a relaxation parameter. Aside from the benefit of an instanteous 
balance in mass flow rate, another merit of this device is that fewer 
global iterations are required in the relaxation solution of equation (15), 
as it is now more nearly satisfied at the outset. 

However, this approach was found def.;ctive, in theory as well as in 
fact. The solving of equation (15) in the form indicated yields nonphysical 
results, as it ;;>rovides a quasi-full potential transonic flow equation 
whose elliptic-hyperbolic transition point can differ markedly from Mach 1. 

This difficulty can be largely alleviated, although not totally circumvented, 
by replacing equation (15) with the equation 

V . (pv^)^ « -7 • (Pu)j^_j (17) 

whose point of transition mere closely approximates the physics of the 
flow. 


Pressure gradients in Dodge's unrevised method would be computed on 
pass n from the equation 





(18) 


In the revised version, pressure gradients are allowed to develop 
during the mass-balancing iteration according to the equations 




3x 

^k,n 

3pn,k 


3y 

^k,n 

3pn.lt 


3z 

-"k,n 



k 

gx 







y yy 




(19) 

( 20 ) 
( 21 ) 


The quantity g is determined through equation (16), and g^ by second- 
order backward differencing. This precedure represents a weak decoupling 
of the streamwise pressure gradient, since the g terms are the dominant 
contributions, and since these contributions are determined from local 
plane-to-plane continuity considerations, somewhat independently of the 
output from the global continuity equation [eq. (17)] on the previous pass. 


7 


NUMERICAL ANALYSIS 


The algorithm deemed most efficient for numerically solving equation 
(17) on the array-processing computer is the Zebra algorithm of South 
et al. (ref. 13). This 3-D relaxation technique is in some respects 
similar to the hopscotch method of Gourlay (ref. 14). In equation (17) 
central differences are applied to all derivative terms. Variables in 
plane i are updated in checkerboard fashion, plane by plane in a 
downstream sweep, using already updateo values at plane i - 1 and old 
iteration values in plane i ■«' 1. Iterative repetition of downstream 
sweeps is used to converge the field, with a relaxation parameter 
employed to speed cdnvergence. 

A second-order accurate, implicit linear multistep method is used 
on equation (11) to march in the stepwise direction. The implicit 
equations are iteratively solved using a checkerboard successive overrelaxa- 
tion scheme, with mass balancing built in as previously described. Stream- 
wise derivatives are backward differenced second-order accurate, while 
derivatives in the (transverse) cross-plane are approximated second-order 
using central differences. A prediction of form 

^i * ^^i-1 ■ ^i-2 (22) 

is used to initially estimate a velocity variable in plane i. The 
checkerboard method is then employed on the differenced counterpart of 
equation (11) to update variables in plane i in two cycles, with values 
updated on cycle i fed into the succeeding cycle. This two-cycle update 
process is iterated, employing equations (16) and (19) to (21) to alter 
the flow speed and pressure gradients until a balance in mass flow is 
achieved. 


DEVELOPING FLOW IN A STRAIGin DUCT 

The revised method of Dodge has been employed to develop a finite- 
difference numerical model for three-dimensional viscous flows in confined 
regions. For boundary- layer resolution, the capability to allow individual 
coordinate stretc.hing in each coordinate direction has been incorporated. 


ORIGINAL PAGE lb 
OP POOR QUALITY 


8 


The metlkod so developed has been programmed using the SLl vector language 
for the Cyber 203 array processor, and appears debugged. The 32-bit 
half-word option of SLl has been employed in programming the Zebra relaxation 
algorithm for solving equation (17), while 64-bit full-word arithmetic is 
used in programming the checkerboard marching algorithm. The program has 
been tested by application to the problem of computing the steady developing 
flow in a straight duct (see fig. 1). Boundary conditions for the problem 
are now given. 

Boundary Conditions 

Inflow: T * Tp - -j u^ 

specified velocity profiles, W ■ H(y,z), 

• R(y»z)» * constant 

^jj(0,y,z) * g(0) • 

Duct walls; velocity no slip, * 0, T * T , P ■ 

II w w 

Outflow: py extrapolated, g^ extrapolated 

Artificial barriers: The computational domain is taken to be one 

quarter of the total duct cross section, and symmetry conditions 
are applied at the two resulting (nonwall) artificial bairiers. 

Here the normal velocity component vanishes together with normal 
derivatives of 4 and the other velocity variables. The variables 
P and T, of course, depend on ♦ and velocity at these boundaries. 
However, for constant total temperature, vanishing normal derivative 
in T, p is the natural bounuary condition. 

On the first pass, the value is allowed to develop in the 
calculation from mass flow rate balancing down the duct. Thereafter, it 
is held fixed. 


COMPUTATIONAL RESULTS 

The method developed and programmed has been exercised by application 
to the duct flow calculation with Reynolds number > 100 and Mach number 
ranging from 0.1 to 0.3 down the duct. An initial pressure distribution 

X 

<ti a ^g(x)dx (23) 

is determined by mass balancing on the first pass. On successive 
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marching passes, this distribution is corrected according to equation (17) 
plus whatever g corrections are necessary in order to balance mass at each 
trossplane of the calculation. The above process is iteratively repeated 
until the maximum change in the pressure field becomes sufficiently small. 

Figures 2 to 4 exhibit computational results after 32 passes. However, 
it can be inferred from the numerical performance indicators shown by 
figures 5 to '7* that the iteration has sufficiently converged to produce 
essentially the same results in around 20 iterations. This figure corres- 
ponds roughly to that given by Dodge (ref. 9) as the number necessary to 
converge a similar problem. However, no real comparison of iteration 
counts to convergence can be drawn, as tolerance levels for convergence of 
various iterative processes could be expected to affect the number of 
passes needed. 

Since the flow was not started at the channel entrance, but with a 
velocity profile supplied at some point farther downstream, not much can 
be said in terms of quantitative assessment of the numerical results, 
which for the most part appear qualitatively excellent. The direction of 
the crossflow and other features in figure 2 appear reasonable and agree 
with that of a computational experiment by Baker (ref. IS) for laminar 
comer flow. The approximately linear variation in centerline pressure 
exhibited by the graph in figure 3 is certainly reasonable for nearly 
incompressible flow, as also exhibited away from th'.; channel entry 
region in a computational experiment of Briley (ref. 1) . Perhaps the 
most questionable feature of the results is the tail-off in centerline 
velocity of figure 4. This could be caused by a calculation not yet 
completely converged near the outflow. However, it is perhaps more 
likely to be the result of the outflow boundary condition treatment. 

For example, the condition 4 ■ constant in the outflow plane forces 4^ 
and 4^ to vanish, which may not be physically reasonable throughout 
the outflow plane, particularly since this does not happen upstream. 


♦Ordinate values in figures 4 to 7 have been magnified by a factor of 

10 . 
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SUM4ARY AND CONCLUSIONS 


A revised version of Dodge's split-velocity method for numerical 

solution of compressible confined flow has been developed. Prelininary 

results for lc‘4 Mach number flow appear encouraging. However, the method 

in general is by no means fully understood or confidently tested. A 

curious feature of the present approach is exhibited by the need for weak 

decoupling of streamwise pressure gradients in order to achieve a 

convergent numerical process. However, Spalding (ref. 12) alleges that, 

in order to achieve physical solutions, a full decoupling is necessary 

with parabolized equations, and Briley (ref. 1) reports successful and 

meaningful calculation obtained using an algorithm incorporating this 

« 

practice. Other questions which bear investigation concern the performance 
of the revised algorithm for higher Reynold's and Mach number flows. To 
gain further ccnfidence in the method, detailed comparisons with independent 
computational results need to be initiated. 

The revision of Dodge's method reported herein is new to the method, 
although classical in physical origins and certainly used previously 
with other computational methods. This investigation has proven the 
checkerboard iteration to be a convenient method for solving implicit 
finite-difference models of the Navier-Stokes equations on the vector 
computer. Further evidence of the computational utility of the zebra 
algorithm for solving the full -potential equation (with a forcing term 
added) in three dimensions has also been gained. 

In conclusion, it is expected that forthcoming investigation will be 
directed: 

(a) to providing details of computation times for the revised Dodge 
method , 

(b) to obtaining comparisons with independent numerical results, and 

(c) to calculating higher Reynolds and Mach number flow. 
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Figure 2. Crossflow and streamwise velocity component at channel exit 


ti 

§e.«| 


6.M 


e^i 


.« 4 14 14 

STREM4ISE DISTANCE 


Figure 3. Centerline pressure variation dovm channel 
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Figure 7 
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Channelwise potential increment Aif vs. iteration count 
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